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The evidence in favor of a dark energy component dominating the Universe, and driving its presently accelerated 
expansion, has progressively grown during the last decade of cosmological observations. If this dark energy is 
given by a dynamic scalar field, it may also have a direct interaction with other matter fields in the Universe, in 
particular with cold dark matter. Such interaction would imprint new features on the cosmological background 
evolution as well as on the growth of cosmic structure, like an additional long-range fifth-force between massive 
particles, or a variation in time of the dark matter particle mass. We review here the implementation of these new 
physical effects in the N-body code GADGET-2, and we discuss the outcomes of a series of high-resolution N-body 
simulations for a selected family of interacting dark energy models, as already presented in Baldi et al. f2Q||. We 
interestingly find, in contrast with previous claims, that the inner overdensity of dark matter halos decreases 
in these models with respect to ACDM, and consistently halo concentrations show a progressive reduction for 
increasing couplings. Furthermore, the coupling induces a bias in the overdensities of cold dark matter and 
baryons that determines a decrease of the halo baryon fraction below its cosmological value. These results go in 
the direction of alleviating tensions between astrophysical observations and the predictions of the ACDM model 
on small scales, thereby opening new room for coupled dark energy models as an alternative to the cosmological 
constant. 



1. INTRODUCTION 

The existence of dark energy (DE) and cold 
dark matter (CDM) as the two most abundant 
components of the Universe keeps acquiring sup- 
port as new and more detailed observational data 
become available. Throughout the last decade, a 
large number of independent and complementary 
datasets (e.g. P, 0, H, 0, have shown that 
the Universe is almost spatially flat, with a total 
matter density of at most 24% of the critical den- 
sity, and the remaining 76% being made of a DE 
component able to drive an accelerated expan- 
sion. The nature of this dominant component, 
however, is still unknown and its understanding 
constitutes a serious theoretical challenge. Dy- 
namical DE models based on the evolution of a 
self-interacting scalar field have been proposed 
d, 0] as a possible way to overcome the severe 
fine-tuning problems that affect the (otherwise 



very successful) cosmological constant. An inter- 
esting idea in the context of dynamical DE models 
is then given by the recently suggested possibility 
of a direct interaction between the dark energy 
scalar field and other matter species in the Uni- 
verse @, i, E [13 ■ Understanding in detail 
what effects this interaction imp rints on observ- 
able features [H [li, [H [H, [H llM bjf is there- 
fore a crucial step for a direct observational test 
of these models. 

In this paper, after briefly recalling the es- 
sential features of interacting DE models and 
discussing their implementation in numerical N- 
body codes, we present the results of the flrst fully 
self-consistent high-resolution hydrodynamic sim- 
ulations of cosmic structure formation for a se- 
lected family of interacting DE cosmologies. 

The results presented here are discussed in 
larg er detail and in a more complete fashion in 
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2. INTERACTING DARK ENERGY 
COSMOLOGIES 

Interacting dark energy cosmologies can be de- 
scribed as a multicomponent system where, al- 
though the conservation of the total stress energy 
tensor T'^i, is not violated: 



0, 



(1) 



the individual stress energy tensor T^^^i^a) of each 
species a is in general not conserved. If so, its 
divergence has a source term Qi^a.)v accounting 
for the coupling between the individual species: 



(2) 



and the conservation of the total stress energy 
tensor (Eqn. [T]) translates into the constraint 



Q(a)i/ = . 



(3) 



We consider then a system described by the La- 
grangian: 

£ = -^d^cj>d^<j> - U{cl)) - m((/.)^7A + rkinM , (4) 

where the role of the DE is played by a scalar 
field (j) with self- interaction potential U{(j>), and 
where the mass of matter fields tp coupled to the 
DE is also a function of (j). The choice to((/)) then 
specifies the source term via the expression: 



91nm(0) 



Pc d^cj) . 



(5) 



Due to the constraint Q, if DE couples only to 
cold dark matter (CDM, hereafter denoted with 
a subscript c) one has that Q(^c)v — ~Q[it>)u- 
the following we will assume that the latter condi- 
tion always holds. In particular, in a cosmological 
system composed by DE, CDM, and baryons, we 
will always assume the baryons (denoted with a 
subscript h) to have no coupling to the DE. 

The zero-component of equation ([2|) provides 
the conservation equations for the energy densi- 
ties of each species: 



Pc 



'3Hp4,{l + w^) 
'"il-Lpb , 



(6) 
(7) 
(8) 



which will determine the background evolution of 
the Universe. Here a prime denotes a derivative 
with respect to the conformal time r, defined via 
dr = dt/a, and Ti. = a' / a is the conformal Hubble 
function, while lu^ = Pcj,/p,f, is the equation of 
state of the DE. 

We will focus our attention here on a class of 
models identified by the choice: 

m(0) = moe-''^(«^^g(,)o = -^^ 



-Pc^', (9) 



where M = l/y/SnGN and Gn is Newton's grav- 
itational constant. 

This set of cosmologies has been widely inves- 
tigated, for the case of a constant coupling func- 
tion, with regard to its background and linear 
perturbation features 0, 0] , to its effects on struc- 
ture formation 21 , 2^ , and also via a first N-body 
simulation [i^ ]. 

By perturbing up to linear order the full system 
of equations governing our cosmological frame- 
work it is possible to derive an effective acceler- 
ation equation for a CDM test particle at a dis- 
tance r from another CDM particle of mass Ale, 
which takes the form (for a complete and rather 
compact derivation see j2( 



Vr. = -H 1- 



/3c (0) 
M H 



■.GrMr. 



(10) 



This equation is manifestly different in three as- 
pects from the usual newtonian acceleration equa- 
tion in an expanding Universe: first, the velocity- 
dependent term now contains an additional con- 
tribution proportional to the DE-CDM coupling 
I3c{4>)\ second, the CDM test particle feels an ef- 
fective gravitational constant Gc given by [IH : 



Gc = Gn [1 + 2/32 (</,)] 



(11) 



third, as a consequence of the energy exchange 
with the DE scalar field (Eqn. [7]), the mass Mc 
of the CDM particle sourcing the gravitational 
potential changes in time according to the expres- 
sion: 



(12) 



In the N-body analysis described in the following, 
we consider (3c{(f>) to be constant, so that the ef- 
fective mass formally reads Mc = Mce"'^=^"*~"*°^ 
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Eqn. [TO] is very important for our discussion 
since it represents the starting point for the im- 
plementation of coupled DE models in an N-body 
code. It is essential to stress here the vectorial 
nature of Eqn. [TOl which is a key point to under- 
stand the effects of the new physics induced by 
the DE-CDM interaction and hence needs to be 
properly taken into account for a correct imple- 
mentation in N-body algorithms. 

3. THE SIMULATIONS 

The main aim of our investigation is to ex- 
plore the effects that arise in the process of cos- 
mic structure formation as a consequence of the 
interaction between the DE scalar field and the 
CDM fluid. Our focus will concentrate on the 
nonlinear regime of structure formation, and to 
this end we study a set of high-resolution cos- 
mological N-body simulations performed with the 
code GADGET-2 [H], that we suitably modified 
for this purpose. Previous attempts to use cosmo- 
logical N-body simulations for different flavours of 
modified Newtonian gravity have been discussed, 
for instance, in 



23| 



is the only previ- 



but to our knowledge 
ous work focusing on the properties of nonlin- 
ear structures in models of coupled quintessence. 
For this reason we apply our methodology to a 
series of interacting DE models identical to the 
ones considered in [23] , where the only difference 
with respect to this previous work amounts to 
updating the cosmological parameters to the lat- 
est WMAP [2] results {flc = 0.213, ni, = 0.044, 
= 0.743, 0-8 = 0.769, h = 0.719, n = 0.963). 
In these models the scalar field cj) has a Ratra- 
Peebles Q self-interaction potential 



A4+0 



(13) 



with a slope a = 0.143, and with a constant 
coupling /3c to CDM particles only, as described 
above; we label them as RPn for values of the 
coupling Pc = n X 0.04, with n ranging from 1 to 
5, in analogy with [2^ 1 . 

For all the models - and for a reference ACDM 
cosmology with the same cosmological parame- 
ters - we run low-resolution simulations in a large 



box (Lbox = 320/1-1 Mpc, Np^.t = 2 x 128^) to 
test the large-scale behaviour of the different cos- 
mologies and the convergence of our implemen- 
tation; for four of these models (ACDM, RPl, 
RP2, RP5) we then run high-resolution simu- 
lations in a smaller box (ibox = 8Qh^^ Mpc, 
^pait = 2 X 512'^) to investigate the properties 
of highly nonlinear structures. 

3.1. Methods 

We now describe one by one the main distinc- 
tive features of our models and their implementa- 
tion in the N-body code GADGET-2. For a more 
detailed description of our numerical approach see 

Modified expansion rate - The interaction between 
DE and CDM modifies the background evolution 
through the presence of an early DE component dur- 
ing all the period of structure formation. The effect 
of such early DE is to modify the expansion history 
of the Universe. For the same set of cosmological pa- 
rameters at 2 = 0, different values of the coupling 
Pc will then lead to different expansion histories that 
will have to be separately computed and then used 
for the N-body time integration. 

Mass variation - As discussed above, the coupled 
species (CDM in our case) feature an effective parti- 
cle mass that changes in time according to Eqn. 1121 
Therefore, the mass of CDM particles within the sim- 
ulation box will have to be accordingly corrected at 
each timestep, while baryon particles will keep a con- 
stant mass. 

Velocity-dependent acceleratton - As shown in 
Eqn. 1101 the coupling induces an additional velocity- 
dependent term in the acceleration equation for CDM 
particles. Therefore, an additional term of the form 



a„ = Pc-rjV 
M 



(14) 



has to be explicitely added to the acceleration of every 
CDM particle at each timestep. 

Fifth-force - Finally, one of the most important 
modifications introduced by the DE-CDM interaction 
is the presence of a long-range fifth-force effectively 
represented by a modified gravitational constant, for- 
mally written as in Eqn. 1111 for the gravitational 
interaction of two CDM particles, while any inter- 
action involving a baryon particle will be governed 
by the standard gravitational constant Gjv. This de- 
pendence of the gravitational strength on the type of 
the particles involved in the interaction requires an 
N-body code to be able to distinguish among difi'er- 
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ent particle types in the gravitational force calcula- 
tion. In GADGET-2, the gravitational interaction is 
computed by means of a TreePM hybrid method (see 
[2^ for details about the TreePM algorithm), so that 
both the tree and the particle-mesh algorithms have 
been modified in order to account for this effectively 
species-dependent gravity. 

3.2. Tests 

As a test of our implementation we check whether 
the linear growth of density fluctuations in the simula- 
tions is in agreement with the linear theory prediction 
for each coupled DE model under investigation. To 
do so, we compute the growth factor from the low- 
resolution simulations and we compare it with the 
solution of the system of equations for the evolution 
of linear perturbations in interacting DE models, nu- 
merically integrated with a suitably modified version 
of the Bohzmann code CMBEASY The compar- 
ison is shown in the left panel of Fig. [1] The accuracy 
of the linear growth computed from the simulations in 
fitting the theoretical prediction is of the same order 
for all the values of the coupling, and the discrepancy 
within the two never exceeds a few percent. 

4. RESULTS 

We present here the most relevant outcomes of the 
high-resolution simulations of coupled DE cosmolo- 
gies described above. 

As first basic analysis steps we apply the Friends- 
otPriends (FoF) and SUBFIND algorithms '34 to 
identify groups and gravitationally bound subgroups 
in each of our simulations. Given that the seed used 
for the random realization of the power spectrum in 
the initial conditions is the same for all the differ- 
ent runs, structures will form roughly at the same 
positions in all simulations. Therefore it is possi- 
ble, in general, to identify the same objects in all the 
simulations and to directly compare their properties. 
We restrict the comparative analysis of our simula- 
tions to the 200 most massive halos identified by the 
FoF algorithm, which have virial masses ranging from 
4.64 X 10^^ h-^MQ to 2.83 x W^^h'^MQ. 

4.1. Halo density profiles 

With a suitable procedure we can select, among 
the 200 most massive halos found at z = in each of 
our simulations, those objects that can be identified 
with certainty as being the same structure emerging 
in the different cosmologies. This leaves 74 objects 
for our comparison analysis. For these halos we com- 



pute the spherically averaged density profiles of CDM 
and baryons as a function of radius around the posi- 
tion of the particle with the minimum gravitational 
potential. 

Interestingly, the halos formed in the coupled DE 
cosmologies show systematically a lower inner over- 
density with respect to ACDM, and this effect grows 
with increasing coupling. This is clearly visible in the 
right panel of Fig. [1] where we show the density pro- 
files of CDM and baryons for the most massive halo 
in our sample in the four cosmologies with different 
couplings for which we have high-resolution simula- 
tions. We remark that this result is clearly incompati- 
ble with the essentially opposite behaviour previously 
reported by [23l |. 

In fact, unlike found in all the halos in our 

comparison sample have density profiles that are well 
fit by the NEW fitting function ^ 



p{r) 



Pcrit 



5* 



(r/r»)(l + r/r,)2 



(15) 



independent of the value of the coupling. The scale 
radius increases for each halo with increasing cou- 
pling l3c, and becomes larger than that found in 
ACDM. Therefore, the coupling does not affect the 
overall shape of the density profiles, but rather the 
transition point between the two asymptotic slopes 
of the NEW function. In other words, the halos be- 
come less concentrated with increasing coupling. 

4.2. Halo concentrations 

For all the 200 most massive halos found in each 
of our high-resolution simulations we compute halo 
concentrations as c = 7-200 / rs , based on our NEW fits 
to the halo density profiles. Here 7-200 is the radius 
enclosing a mean overdensity 200 times the critical 
density. 

Consistently with the trend found for the inner 
overdensity in the halo density profiles, we find that 
halo concentrations are on average significantly lower 
for coupled DE models with respect to ACDM, and 
the effect again increases with increasing coupling Pc- 
This behaviour is shown explicitly in the left panel of 
Fig. O where we plot halo concentrations as a func- 
tion of the halo virial mass M200 for our four high- 
resolution simulations. 

It is possible to show (see for details) that this 
effect is not due to a change in the formation epoch 
of structures in the coupled cosmologies with respect 
to ACDM. Indeed, a more detailed investigation of 
this effect shows that the decrease of concentrations 
arises as a consequence of the expansion of virialized 
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Fi gure 1 . Left Panel: Evolution of the growth function with redshift for the models of coupled DE investigated with the 
low-resolution simulations. The solid lines are the total growth functions as evaluated numerically with CMBEASY, while 
the triangles are the growth function evaluated from the simulations. Right Panel: Density profiles of CDM (solid lines) 
and baryons (dot-dashed lines) for the most massive halo in the simulation box at 2 = 0. The vertical dot-dashed line 
indicates the location of the virial radius for the ACDM halo. 



structures. This expansion is a direct consequence 
of the DE-CDM interaction, and can be explained 
as follows: on one side, the potential wells of CDM 
halos become shallower as time goes by as a conse- 
quence of the mass decrease of their CDM content, 
as described by Eqn. 1121 on the other side, the new 
velocity-dependent acceleration (Eqn. I14p speeds up 
the motion of CDM particles, thereby injecting en- 
ergy into the system. The combination of these two 
effects determines an overall increase of the total en- 
ergy of the system, which will then slightly expand to 
restore virial equilibrium. 

It is interesting to note that the effects we find go 
in the direction of less "cuspyness" of halo density 
profiles, which is preferred by observations and thus 
in fact opens up new room for the phenomenology of 
interacting DE models. 

4.3. Halo baryon fraction 

The extra force proportional to (see Eqn. Ilip 
that affects only the interaction between two CDM 
particles induces a bias in the evolution of density 
fluctuations of baryons and CDM [H, IH, IH . This 
bias will appear at all scales during the period of lin- 
ear growth of density fluctuations, which makes it 
clearly distinguishable from the hydrodynamic bias 
arising only at small scales as structure evolves. Fur- 
thermore, as it was shown in 23] and 20], this grav- 



itational bias is enhanced in the nonlinear stages of 
gravitational clustering, and reaches its largest am- 
plitude in the core of highly overdense regions. 

It is interesting that the above effect produces a 
baryon deficit in virialized halos, i.e. they contain 
fewer baryons than expected based on their mass and 
the universal cosmological baryon fraction. In partic- 
ular, this means that one can not expect that baryon 
fractions determined through X-ray measurements in 
clusters would yield the cosmological value. In or- 
der to give a rough estimate of the magnitude of the 
discrepancy we compute the relative baryon fraction 
within the virial radius r2oo of all the 200 halos in our 
sample, defined as 



Yb 



where 



A4(<'-2oo) 

Aftot(<r200) 



(16) 



For the ACDM case, our results for the evolution of Yi, 
are consistent with the value of Yi, ~ 0.92 found by the 
Santa Barbara Cluster Comparison Project [s^, and 
with the more recent results of [13] and [4j , while for 
the coupled models the relative baryon fraction shows 
a progressive decrease with increasing coupling, down 
to a value of Yj, ~ 0.86 — 0.87 for the RP5 case, as 
shown in the right panel of Fig. [2] 

It is also important to notice that this effect is al- 
ways towards lower baryon fractions in clusters with 
respect to the cosmological value. This could in fact 
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Figure 2. Left Panel: Evolution of the mean halo concentration as a function of mass for the 200 most massive halos in 
our simulations and for the different cosmological models under investigation. The halos have been binned by mass, and 
the mean concentration in each bin is plotted as a filled circle. Right Panel: Evolution with virial mass M200 of the relative 
baryon fraction Yi, within the virial radius r200 for all the halos in our sample. The coloured diamonds represent the relative 
baryon fraction of each single halo, while the filled circles and the coloured curves show the behaviour of the mean relative 
baryon fraction in each mass bin for our four high-resolution simulations. 



alleviate tensions between the high baryon abundance 
estimated from CMB observations, and the somewhat 
lower values inferred from detailed X-ray observations 
of galaxy clusters [H, [H, H . 

5. CONCLUSIONS 

We have investigated the effects that arise in the 
nonlinear regime of structure formation in the con- 
text of interacting DE models, by means of detailed 
high-resolution N-body simulations run with a suit- 
ably modified version of the code GADGET-2. The 
numerical implementation we have developed is quite 
general and not restricted to the simple specific mod- 
els of coupled quintessence that we have discussed in 
this work. Instead, it should be well suited for a much 
wider range of DE models. 

We have presented here the main conclusions of our 
analysis, concerning the properties of collapsed struc- 
tures in interacting DE models. In particular, we have 
shown that CDM halo density profiles are remarkably 
well fit over the resolved range by the NEW formula 
for any value of the coupling, but there is a clear 
trend of a decrease of the inner halo overdensity with 
respect to ACDM with increasing coupling (or, equiv- 
alently, an increase of the scale radius for increasing 
coupling). This result conflicts with previous claims 



for the same class of coupled DE models [i^]. Con- 
sistently, also halo concentrations are reduced with 
increasing coupling with respect to ACDM. 

Finally, as already shown in we also find that, 
as a consequence of the different effective gravita- 
tional strength between CDM and baryons, these two 
components develop a bias in the amplitude of their 
relative density fluctuations, and this bias is enhanced 
in the nonlinear region within and around massive ha- 
los. The enhancement of the bias in highly nonlinear 
structures has an impact on the determination of the 
baryon fraction from cluster measurements, and we 
have computed for all our halos its evolution with 
coupling, finding that the baryon fraction is reduced 
with increasing coupling by up to ~ 8 — 10% with 
respect to ACDM for the largest coupling value we 
consider. 
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